﻿#pragma once

#include <xy/geom/manifold.h>

namespace xy
{

namespace geom
{

template <std::size_t N> class pbspline : public manifold<2, N>
{
    using vecd_t = vec<2, double>; // 参数域向量
    using vec_t = vec<N, double>;  // P

  protected:
    std::vector<vec_t> m_cpts;
    std::array<std::vector<std::vector<double>>, 2> m_knots;

  public:
    /**
     * @brief pbspline 构造函数
     *
     * @param[in] cpts
     * @param[in] knots
     */
    pbspline(std::vector<vec_t> &&cpts, std::array<std::vector<std::vector<double>>, 2> &&knots)
        : m_cpts(std::move(cpts)), m_knots(std::move(knots))
    {
        if (m_cpts.size() != knots[0].size() || m_cpts.size() != knots[1].size())
            throw std::logic_error("Control points not match.");
    }

    /**
     * @brief 获得控制点数
     *
     * @return std::size_t
     */
    std::size_t pole_num() const
    {
        return m_cpts.size();
    }

    /**
     * @brief 获得控制点
     *
     * @param[in] i
     * @return vec_t
     */
    vec_t pole(std::size_t i) const
    {
        return m_cpts[i];
    }

    /**
     * @brief 设置控制点
     *
     * @param[in] p
     * @param[in] i
     */
    void set_pole(const vec_t &p, std::size_t i)
    {
        m_cpts[i] = p;
    }

    /**
     * @brief 获得控制点集
     *
     * @return decltype(m_cpts)
     */
    auto poles() const -> const auto &
    {
        return m_cpts;
    }

    /**
     * @brief 移动控制点集
     *
     * @return decltype(m_cpts)
     */
    auto move_poles() -> decltype(m_cpts)
    {
        return std::move(m_cpts);
    }

    /**
     * @brief 设置控制点集
     *
     * @param[in] cpts
     */
    void set_poles(std::vector<vec_t> &&cpts)
    {
        m_cpts = std::move(cpts);
    }

    /**
     * @brief 获得 i 维度上的次数
     *
     * @param[in] i
     * @return std::size_t
     */
    std::size_t degree(std::size_t i = 0) const
    {
        return m_knots[i].size() - 2;
    }

    /**
     * @brief 获得每个维度上的次数
     *
     * @return std::array<std::size_t, 2>
     */
    const std::array<std::size_t, 2> &degrees() const
    {
        return {degree(0), degree(1)};
    }

    /**
     * @brief 获得 i 维度上的节点数
     *
     * @param[in] i
     * @return std::size_t
     */
    std::size_t knot_num(std::size_t i = 0) const
    {
        return m_knots[i].size();
    }

    /**
     * @brief 获得 i 维度上的节点向量
     *
     * @param[in] i
     * @return const std::vector<std::vector<double>>&
     */
    const std::vector<std::vector<double>> &knots(std::size_t i = 0) const
    {
        return m_knots[i];
    }

    /**
     * @brief 设置 i 维度上的节点向量
     *
     * @param[in] i
     * @return std::vector<std::vector<double>>&
     */
    std::vector<std::vector<double>> &knots(std::size_t i = 0)
    {
        return m_knots[i];
    }

    /**
     * @brief 移动 i 维度上的节点向量
     *
     * @param[in] i
     * @return std::vector<std::vector<double>>
     */
    std::vector<std::vector<double>> move_knots(std::size_t i = 0)
    {
        return std::move(m_knots[i]);
    }

    /**
     * @brief 计算曲线上点的位置
     *
     * @param[in] param
     * @return vec_t
     */
    vec_t d0(const vecd_t &param) const override
    {
        vec_t ret;
        for (std::size_t i = 0; i < m_cpts.size(); i++)
        {
            double coeff = 1.0;
            for (std::size_t j = 0; j < 2; j++)
                coeff *= one_basis_func(m_knots[j][i].data(), m_knots[j][i].size(), degree(j), 0, param[j]);
            ret += m_cpts[i] * coeff;
        }
        return ret;
    }

    /**
     * @brief 计算曲线上点的位置导数
     *
     * @param[in] param
     * @return std::array<vec_t, 3>
     */
    std::array<vec_t, 3> d1(const vecd_t &param) const override
    {
        std::array<vec_t, 3> ret;
        for (std::size_t i = 0; i < m_cpts.size(); i++)
        {
            auto uders = deriv_one_basis_func(m_knots[0][i].data(), m_knots[0][i].size(), degree(0), 0, 1, param[0]);
            auto vders = deriv_one_basis_func(m_knots[1][i].data(), m_knots[1][i].size(), degree(1), 0, 1, param[1]);
            ret[0] += m_cpts[i] * uders[0] * vders[0];
            ret[1] += m_cpts[i] * uders[1] * vders[0];
            ret[2] += m_cpts[i] * uders[0] * vders[1];
        }
        return ret;
    }

    /**
     * @brief 计算曲线上点的位置二阶导数
     *
     * @param[in] param
     * @return std::array<vec_t, 6>
     */
    std::array<vec_t, 6> d2(const vecd_t &param) const override
    {
        std::array<vec_t, 6> ret;
        for (std::size_t i = 0; i < m_cpts.size(); i++)
        {
            auto uders = deriv_one_basis_func(m_knots[0][i].data(), m_knots[0][i].size(), degree(0), 0, 2, param[0]);
            auto vders = deriv_one_basis_func(m_knots[1][i].data(), m_knots[1][i].size(), degree(1), 0, 2, param[1]);
            ret[0] += m_cpts[i] * uders[0] * vders[0];
            ret[1] += m_cpts[i] * uders[1] * vders[0];
            ret[2] += m_cpts[i] * uders[0] * vders[1];
            ret[3] += m_cpts[i] * uders[2] * vders[0];
            ret[4] += m_cpts[i] * uders[0] * vders[2];
            ret[5] += m_cpts[i] * uders[1] * vders[1];
        }
        return ret;
    }

    /**
     * @brief 获得参数范围
     *
     * @param[in] i
     * @return std::array<double, 2>
     */
    std::array<double, 2> bound(std::size_t i = 0) const override
    {
        auto min = std::numeric_limits<double>::max();
        auto max = -std::numeric_limits<double>::max();
        for (const auto &knot : m_knots[i])
        {
            min = std::min(min, knot.front());
            max = std::max(max, knot.back());
        }
        return {min, max};
    }
};

} // namespace geom

} // namespace xy